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Abstract.  A  recent  multigrid  method  for  computing  steady  inviscid  compressible  flow  is  analyzed 
for  the  one-dimensional  scalar  case.  The  discretization  in  space  by  means  of  upwind  differencing 
has  first-  or  second-order  accuracy.  Three  types  of  relaxation  schemes  that  use  only  local  informa¬ 
tion  are  examined.  A  numerical  experiment  with  a  smooth  steady  solution  shows  good  agreement 
with  the  estimated  damping  rates.  A  discontinuous  solution  displays  slower  convergence.  This  is 
mainly  caused  by  the  singularity  at  the  shock.  The  singularity  can  be  removed  by  the  additional 
constraint  of  local  conservation,  which  is  achieved  through  a  special  kind  of  local  relaxation,  using 
information  about  conservation  from  coarser  grids.  The  sonic  point  causes  some  slow-down  as 
well,  but  this  can  be  neglected  in  practical  applications.  It  turns  out  that  a  first-order-accurate 
solution  can  be  obtained  by  1  F-cycle  per  grid,  whereas  second-order  accuracy  requires  about  4. 

Key  words:  multigrid  method,  hyperbolic  conservation  laws,  shock  waves 
1.  Introduction 

The  multigrid  technique  is  a  powerful  tool  for  the  construction  of  O(N)  methods  for  a  wide 
class  of  problems.  Especially  the  solution  of  elliptic  PDE’s  has  received  much  attention,  and 
the  technique  is  well  established  for  these  problems  [3,4,13].  It  has  taken  some  time  before  any 
results  for  hyperbolic  equations,  specifically  the  Euler  equations  of  gas-dynamics,  were  obtained, 
but  successful  experiments  are  now  available  [6,7,10]. 

In  [10]  it  was  demonstrated  that  the  O(N)  behaviour  can  indeed  be  obtained  for  the  relatively 
simple  problem  of  transonic  flow  through  a  channel  with  a  circular  bump  on  one  wall.  The 
multigrid  technique  used  is  essentially  standard,  thanks  to  discretization  by  means  of  upwind 
differencing.  It  is  applied  to  a  linear  system  based  on  the  nonlinear  equations.  A  direct  nonlinear 
formulation,  known  as  a  Full  Approximation  Storage  (FAS)  scheme,  is  not  fundamentally  different 
from  this  approach.  Indeed  it  was  shown  in  [5]  that  the  nonlinear  version  of  the  technique  in  [10], 
with  minor  variations,  yields  a  similarly  good  performance.  FAS  is  more  expensive  [7],  but  more 
versatile,  especially  when  local  grid-refinement  is  used. 

A  major  drawback  of  the  technique  so  far,  is  the  necessity  to  use  many  multigrid  cycles,  as 
opposed  to  the  one,  or  few,  cycles  that  are  usually  sufficient  for  elliptic  problems.  An  example  can 
be  found  in  [11].  There,  however,  an  ill-posed  problem  was  studied,  and  the  solution  was  forced 
to  convergence  far  below  the  truncation  error  (which  could  not  be  completely  defined  because  of 
the  ill-posedness).  But  even  under  better  circumstances  the  multigrid  technique  does  not  seem 
to  completely  live  up  to  the  expectations  [6,7,8]. 

This  paper  provides  some  insight  in  what  can  actually  be  expected,  and  what  is  the  cause  of  - 
failing  expectations,  by  considering  the  simplest  nonlinear  model  equation  one  can  think  of:  the 
one-dimensional  scalar  inviscid  Burgers’  equation. 


This  work  ha*  been  supported  by  the  Center  for  Large  Scale  Scientific  Computing  (CLaSSiC)  Project  at 
Stanford  under  the  Office  of  Naval  Research  Contract  N00014-82-K-0335. 
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Both  a  smooth  and  a  discontinuous  test  problem  are  considered  ($2).  The  spatial  discretiza- 
tion  is  obtained  through  upwind  differencing  by  Roe’s  scheme  [12],  with  first-  or  second-order 
accuracy. 

In  §3  three  relaxation  schemes  are  described  that  use  only  local  information.  This  property 
is  desirable  for  vectorization,  parallelization,  and  local  grid- refinement.  Their  smoothing  rates 
are  estimated  from  a  linearized  version  of  the  nonlinear  differential  equation. 

The  multigrid  scheme  is  described  in  §4.  By  combining  the  relaxation  operator  (§3.2)  with 
the  Coarse-Grid  Correction  (CGC)  operator  (§4.1),  we  can  obtain  the  asymptotic  convergence 
rates  for  various  schemes  (§4.2).  Some  ideal  relaxation  schemes,  that  have  an  estimated  damping 
rate  of  zero  in  regions  where  the  velocity  is  well  above,  or  well  below,  zero,  are  considered  in  §4.3. 
The  Full  Multigrid  scheme  (FMG),  with  its  restriction  and  prolongation  operators,  is  described 
in  §4.4. 

The  associated  amount  of  work  is  estimated  in  §4.5.  Solution  errors  and  stopping  criteria 
are  reviewed  in  §4.6.  The  practical  construction  of  a  relaxation  scheme  is  described  in  §4.7. 

The  analysis  breaks  down  if  the  velocity  of  the  flow  changes  its  sign.  At  the  shock,  the  discrete 
equations  become  singular.  This  is  appropriate,  because  the  original  differential  equations  are 
singular.  The  singularity  can  be  removed  by  augmenting  the  differential  equations  with  the 
Rankine-Hugoniot  jump  relations,  which  are  implicit  in  the  conservation  form  of  the  equations. 
The  time-accurate  discrete  equations  in  conservation  form  resolve  the  singularity  in  the  same  way. 
However,  time-accuracy  has  been  sacrificed  here  to  obtain  fast  convergence.  Thus,  we  need  an 
alternative  way  to  remove  the  singularity.  This  can  be  achieved  by  imposing  local  conservation 
through  a  special  kind  of  local  relaxation,  using  information  about  conservation  from  the  coarser 
grids  (§4.8). 

At  the  sonic  point,  we  analyze  the  damping  rate  in  a  different  way,  now  using  the  linearized 
difference  equations.  Only  Point-Jacobi  for  a  first-order-accurate  solution  is  considered  here 
(§4.9). 

Numerical  Tesults  on  the  two  test  cases  are  presented  in  §5.  A  comparison  between  estimated 
and  experimental  asymptotic  damping  rates  is  made  (§5.1).  Next,  the  minimum  number  of  cycles 
required  to  obtain  a  steady  state  within  the  truncation  error  is  considered  (§5.2).  Two  multigrid 
extensions,  the  deferred  correction  technique  and  r-extrapolation,  are  discussed  in  §5.3  and  §5.4, 
respectively. 

The  main  conclusions  are  summarized  in  §6. 

2.  Test  Problems 

We  consider  two  types  of  steady  solutions  to  the  one-dimensional  scalar  hyperbolic  equation 

««  +  (£«*)»  =*(*),  (2.1) 

the  first  being  smooth,  the  second  discontinuous.  These  solutions  are  periodic  in  x  on  the  interval 
[0, 1).  Setting 

u(x)  =  Co  +  sin2ir(*  —  (),  Co  =  2,  (2.2a) 

leads  to  a  source  term 

*(*)  =  2r  (c o  +  sin  2x(x  —  {))  cos  2x(x  —  ().  (2.26) 

Although  the  steady  state  corresponding  to  this  source  term  is  continuous,  a  shock  may  occur 
during  the  evolution  toward  the  steady  state. 

A  discontinuous  steady  solution  is  found  for 


s(x)  =  sin  2t(z  -  £), 


(2.3a) 
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namely, 

«(*)  = 

For  0  <  £  <  this  solution  has  a  sonic  point  at  *  =  (  and  a  shock  at  x  =  £  +  5  (cf.  [15]). 
Numerical  examples  will  be  given  for  £  =  0.1.  In  both  cases  the  solution  is  periodic.  It  should  be 
remarked  that  Eq.(2.1)  must  be  interpreted  as  the  limit  for  zero  viscosity  v  of  the  same  equation 
with  vuxx  added  to  the  right-hand  side. 

A  computational  grid  with  N  zones  is  defined  by  *t  —  (k  +  |)/i,  where  the  cell-size  h  =  l/N 
and  k  —  0, . . . ,  N  —  1.  Equation  (2.1)  is  discretized  in  space  by  averaging  per  volume: 

=Ihu  =  l  J  dx  u(x),  (2.4) 

The  same  discretization  operator  lh  is  applied  to  the  residual 

r(tt,z)  =  «(r)-(luJ)„  (2.5) 

yielding 

ri(u)  =  *i  -  ^(/t+i  (2-6) 

The  fluxes  /t±  ^  are  evaluated  by  upwind  differencing .  Godunov’s  scheme  lets 

/*+$  =/(«  *,«*+i)  =  £max(max(0,ut)J,min(0,ut+i)J).  (2.7) 

We  use  Roe’s  scheme  [12,15],  which  in  this  case  is  identical  to  Godunov’s  scheme  except  for  its 
treatment  of  the  sonic  point: 

fk+i  =/(«*. «*+i)=  if  «i  <  0  <  ut+i.  (2.8) 

With  this  flux,  the  solution  at  the  sonic  point  is  smooth,  whereas  Godunov’s  scheme  sets  /t+^  =  0, 
causing  an  O(h)  jump. 

The  numerical  steady  solution  u*  obeys  rt(ut_i, u*,!<t+1)  =  0  for  »  =  0, ...,1V  —  1.  It  is 
first-order  accurate,  i.e.,  ||uj  —  ujt||  =  0(h),  where  u*t  =  lhu(x)  is  the  average  per  cell  of  the 
exact  solution,  and  ||  •  ||  is  a  relevant  norm. 

Second-order  spatial  accuracy  can  be  obtained  by  assuming  the  solution  to  be  piecewise 
linear  rather  than  piecewise  constant  [15].  The  idea  is  to  write  the  solution  as 

«(*)  =  „* +(!=.£*)**,  !!=-£*!<!.  (2.9) 

Here  At  is  related  to  the  average  derivative  of  the  solution  within  cell  k: 

ik+^h 

At  S  J  dx  g.  (2.10) 


{ 


sinx(x  —  £) 
—  sin  r(x  —  £) 


for  0  <  x  <  £  +  5, 
for  £  +  i  <  x  <  1. 


(2.36) 


Its  numerical  value  can  be  computed  from  ut-i,ut,ut+i  by 


At  =  ave(uk  -  ut_i,  «t+i  -  ut). 


(2.11) 
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In  smooth  regions  we  want  at>e(A_ ,  A+)  =  |(A_  +  A+),  whereas  some  limiting  to  the  smaller  of 
A-  and  A+  is  required  near  discontinuities  to  preserve  monotonocity  (i.e.,  to  avoid  local  under- 
or  overshoots).  We  use  an  averaging-limiting  procedure  from  [1]: 


/A  A  x  (Ai  +  c2)A-  +  (A2_+e2 

“”<4-  •  4*> - ai  +  4+2  <1 


The  bias  ca  prevents  division  by  zero.  Clipping  of  smooth  extrema  is  avoided  if  c„  approxi¬ 
mately  equals  the  average  value  of  |u*  —  ut_i|  in  smooth  regions  of  the  flow.  For  the  numerical 
examples  mentioned  above  we  adopt  ea  =  4  A  for  the  smooth  problem  (2.2)  and  ca  =  2h  for 
the  discontinuous  problem  (2.3).  Once  the  A*  are  computed,  the  fluxes  /i+^  are  evaluated  by 
/(tut  +  iAt,  tit+i  -  iAi+i),  using  either  Godunov’s  (2.7)  scheme  or  Roe’s  scheme  (2.7, 2.8).  It 
should  be  remarked  that  the  scheme  has  second-order  accuracy  almost  everywhere.  In  the  worst 
case,  there  remains  an  0(h)  error  in  the  cells  near  the  shock.  The  rest  of  the  solution  has  an 
0(h 2)  error. 

Apart  from  the  smoother  behaviour  of  the  first-order  Roe  scheme  near  the  sonic  point,  there 
is  another  reason  for  preferring  Roe’s  scheme  above  Godunov’s,  as  will  be  seen  in  §4.7. 

The  numerical  problem  of  determining  the  steady  state  to  (2.1)  can  be  written  as 

ri(ut_p,...,ut+p)  =  0,  for  fc  =  0,  (2-13) 

Here  p  =  1  for  a  first-order  and  p  =  2  for  a  second-order  scheme.  For  the  example  given  in 
Eq.(2.3),  (2.13)  is  ill-posed.  Suppose  that  ul  >  0  is  the  state  in  the  cell  at  the  left  from  the  shock, 
ur  <  0  at  the  right,  and  that  uj u  is  the  state  in  the  cell  that  contains  the  shock:  ul  >  >  ur. 

Then  the  residual  for  the  first-order  scheme  ru  —  *m  —  (jUr  —  juJ),  which  is  independent  of 
UAr-  Thus,  the  state  u m  never  enters  the  discrete  equations  (2.13)  and  can  be  chosen  arbitrarily, 
as  long  as  ux,  >  UM  >  ur.  For  a  second-order  scheme,  ul  and  ur  are  replaced  by  u l  +  |  A l  and 
ur  —  iA/i  in  the  expression  for  the  residual.  The  average  differences  Ajt  and  A r  will  depend 
on  u m,  but  if  the  latter  state  is  not  close  to  either  ul  or  ur,  limiting  will  occur.  In  this  way 
the  A l  and  A r  become  small,  and  the  residual  behaves  practically  the  same  as  in  the  first-order 
case.  As  a  consequence,  one  expects  the  numerical  non-uniqueness  to  persist  in  the  second-order 
scheme,  although  um  will  be  confined  to  a  smaller  interval  than  that  for  the  first-order  case. 

The  singular  behaviour  at  the  shock  is  the  immediate  consequence  of  dropping  the  time- 
dependency  from  the  original  equation.  In  a  time-accurate  integration,  the  conservation  form  of 
the  equation  automatically  takes  care  of  the  jump  equation  across  the  discontinuity.  Then 


1  N~l 

—  ^  u *  =  C  (a  constant), 


for  each  time-step.  In  this  way  the  singularity  does  not  occur.  If  the  integration  is  not  time- 
accurate,  the  final  tZ*f  can  be  recovered  from 

N- 1 

um  =  NC-  ^  uk,  (2.15) 

isO,*/Af 

but  only  if  there  is  one  shock.  Alternatively,  the  detailed  structure  in  zone  k  =  M  may  be  used 
to  find  u m  (see  [15]).  In  this  paper  we  will  use  local  conservation  as  an  additional  constraint 
to  remove  the  singularity.  Conservation  is  imposed  on  the  coarsest  grid.  On  fine  grids  the 
information  about  conservation  is  obtained  from  the  next  coarser  grid  (§4.8). 


J*  /  J-  -• 


3.  Relaxation  schemes 

3.1.  Construction.  The  problem  of  determining  a  steady  solution  to  Eq/  '3)  can  be  written 
in  a  more  abstract  form  by  introducing  a  nonlinear  operator 

I(u)  =  (^uI),  (3.1a) 

and  its  discrete  version 

Lhk(*k)  =  -/*-*)•  (3-16) 

The  error  vh  =  t<*  —  u\  where  IT*  is  the  steady  state,  obeys 

/,*(„*  +  „»)_£>*)  =  r\  (3.2) 

A  relaxation  scheme  provides  an  approximate  solution  vh  to  Eq.(3.2).  An  example  is  the  explicit 
scheme 


or  the  implicit  scheme 

Is  "*(£)*!  (3-4) 

Here  a  determines  the  amount  of  under-  (a  >  1)  or  over-relaxation  (a  <  1).  The  choice  a  =  1 
yields  a  "backward  Euler”  scheme.  If  the  time-step  At  becomes  large,  the  latter  reduces  to 
Newton’s  method.  This  property  is  exploited  by  the  Switch  Evolution /Relaxation  (SER)  scheme, 
proposed  in  [9,16],  where  At  a  l/||rh||.  After  some  initial  time-dependent  searching  while  the 
residuals  are  large,  this  scheme  automatically  switches  to  Newton’s  method  when  the  solution 
approaches  the  steady  state. 

For  a  one-dimensional  problem,  Eq.(3.4)  with  a  residual  (2.13)  of  first-order  accuracy  (p—  1) 
requires  the  solution  of  a  tridiagonal  system.  This  system  can  be  approximated  by  retaining  only 
the  unaltered  main-diagonal.  Alternatively,  it  can  be  replaced  a  diagonal  system  composed  of  the 
original  main-  and  off-diagonals  [16]. 

For  the  purpose  of  analysis  we  consider  a  linearized  version  of  the  original  equation: 

«t  +  au*  =  *(*)•  (3.5) 

Here  a  is  assumed  to  be  positive,  for  the  moment.  Equation  (3.2)  reduces  to 

Lhvh  =  r\  (3.6a) 

and  a  relaxation  scheme  provides  an  approximate  solution  vh  through 

LhAtvh  =  Lhvh,  A,v  =  vh-vh.  (3.66) 

Here  it  is  assumed  that  the  approximation  Lh  to  Lh  makes  (3.6b)  easy  to  solve.  The  last  equation 
can  be  rewritten  in  terms  of  a  relaxation  operator 


Sh  =  I-(~Lky1Lh 


(37  a) 


Let  the  usual  CFL-number  be  denoted  by  o  =  At  \a\/h  and  the  shift  operators  by  Tut  =  «*+ 1, 
T~lut  =  tii-i-  For  positive  a,  we  have 


Lk  =  £(1  —  T~1). 

A 

(3.8) 

The  explicit  scheme  is  given  by 

*-=• 

(3.9) 

whereas  the  implicit  scheme  lets 

ik  = 

[=-(£)*] 

(3.10) 

a  lower  bidiagonal  system.  Replacing  the  translation  operator  by  a  scalar  multiplication  ui  results 
in 

*fc  =  i?[1  +  a<r(1~w)]’  (3-11) 

a  diagonal  matrix.  For  the  scalar  one-dimensional  problem  studied  here,  the  last  expression 
results  in  the  same  relaxation  operator  as  the  one  for  an  explicit  scheme: 


v-1 

1 

f-. 

1 

f— < 

1 

rd 

II 

■e 

(3.12) 

where 

0  =  <r, 

(3.13a) 

for  an  explicit  scheme,  and 

a_  * 

l  +  o<r(l-w) 

(3.136) 

for  the  diagonal  approximation  (3.11)  of  the  implicit  scheme  (3.4).  The  relaxation  scheme  (3.12) 
can  be  analyzed  for  different  values  of  the  parameter  0.  Although  the  explicit  and  diagonal 
approximation  of  the  implicit  scheme  are  equivalent,  according  to  (3.12),  their  difference  lies  in 
the  construction  of  a  fixed  value  for  0  (3.13).  In  the  case  of  non-constant  a,  or  for  systems  of 
equations,  the  implicit  scheme  provides  more  clues  for  this  construction  then  the  explicit  one.  We 
will  leave  this  subject  now,  and  assume  that  0  is  merely  an  adjustable  parameter  of  the  relaxation 
scheme. 

3.2.  Smoothing  rates.  To  analyze  the  stability  and  damping  of  the  relaxation  schemes,  we 
write  vh  as  a  Fourier-series 

N-l  .  . 

«£  =  £>?exp(2 **-),  k  =  0,-  ■  ,N  —  1.  (3.14) 

i=0 

Let  the  Fourier  operator  Fk  and  its  inverse  (Fk)~1  have  elements 

Fi\  =  ^exp(— 2ri^),  (Fh)~}  =  exp  (2*.^),  (3.15) 

then 

vk  =  Fkvh,  vk  =  (Fh)'1  vk.  (3.16) 

The  relaxation  operator  in  the  Fourier  domain  becomes 

Sk  =  FkSk  ( Fky 


(3.17) 


The  scheme  is  stable  if  the  spectral  radius  p(Sh )  <  1.  The  symbol  of  the  translation  operator  T 
is 

T  =  exp(*0),  6  =  2vl/N.  (3.18) 

For  practical  purposes  6  may  be  assumed  to  be  a  continuous  quantity  (N  — ►  oo).  The  damping 
of  the  scheme  can  be  analyzed  by  evaluating  the  amplification  factor  p(6)  =  Sh(6). 

A  multigrid  method  requires  a  relaxation  scheme  that  has  good  smoothing  for  short  waves. 
We  will  now  consider  three  relaxation  schemes  based  on  (3.12),  which  can  be  easily  generalized 
to  systems  in  one  or  more  dimensions,  use  only  local  information,  and  can  be  vectorized. 

The  first  relaxation  scheme  to  be  considered  is  Point-Jacobi  relaxation  (PJ).  All  cells  in  the 
computational  domain  are  updated  independently.  The  corresponding  growth-factor  is  given  by 

p(0)  =  1  -  0(1  -  f"1).  (3.19) 

Stability  requires 

|/x(0)|2  =  1  -2/3(1  -0)(l-  cos  0)<1,  6  £  [-*,*],  (3.20a) 

implying 

0  <  0  <  1.  (3.206) 

Choosing  0  —  0  does  not  make  sense,  of  course. 

For  a  single-grid  scheme,  the  choice  0=1  has  some  special  properties.  The  growth-factor 
is  n(8)  —  T~i  in  that  case.  The  explicit  scheme  runs  at  the  maximum  local  time-step.  For  an 
implicit  scheme,  a  step  with  0=1  can  be  interpreted  as  a  Point-Jacobi  iteration  for  Newton’s 
method  [PJ  implies  u  =  0,  i.e.,  the  off-diagonals  are  neglected;  Newton’s  method  is  obtained  for 
a  =  1  and  a  very  large  time-step  At,  i.e.,  <r  — ►  oo).  In  a  physical  interpretation,  the  local  errors  can 
be  viewed  as  waves.  For  0  =  1  these  we  convected  over  a  distance  h,  the  cell-size.  For  a  problem 
with  boundaries,  or  with  shocks,  we  thus  expect  a  single-grid  scheme  with  0  =  1  to  yield  fastest 
convergence.  Errors  (in  the  sense  of  deviations  from  the  numerical  steady  state)  are  convected 
as  fast  as  possible  toward  the  boundaries,  where  they  leave  the  computational  domain,  or  toward 
a  shock,  where  they  are  absorbed  by  the  dissipation  in  the  shock.  Because  |^(0)|  s  1,  there 
is  no  damping  in  smooth  regions  of  the  flow,  which  is  consistent  with  the  differential  equation. 
Conservation  in  time  and  time-accuracy,  however,  are  lost.  The  latter  is  not  very  important  for 
steady-state  computations,  as  long  as  the  correct  final  steady  state  is  found.  (For  complicated 
problems,  a  SER  scheme  may  be  safer  [16].)  If  conservation  in  time  is  necessary  for  the  uniqueness 
of  the  final  solution,  it  can  be  imposed  as  a  global  condition  at  the  end  of  one  or  more  relaxation 
sweeps.  For  a  periodic  shockless  problem,  as  in  (2),  the  differential  equation  has  no  mechanism 
to  get  rid  of  the  initial  errors,  and  a  steady  state  can  formally  never  be  reached.  In  practice,  the 
discretization  error  in  the  numerical  flux  provides  a  small  amount  of  dissipation  that  damps  the 
error.  The  relaxation  scheme  can  provide  additional  dissipation.  For  PJ,  the  maximum  amount 
of  dissipation  is  obtained  if  0  =  i. 

The  second  relaxation  scheme  is  a  Multi-Stage  scheme  (MS).  Here  we  consider  only  2  stages. 
The  first  stage  is  a  PJ  step  with  0  =  0i,  advancing  the  solution  u  to  u'.  The  second  stage  is 
again  PJ  relaxation,  but  with  0  =  0?.  Now  the  solution  u  is  advanced  to  u,  using  the  residual 
r'  =  r(u').  The  growth-factor  for  one  MS  step  is 

m  =  1  -  02(1  ~  f-')[l  -  01{1  -  f-1)],  (3.21a) 

which  is  stable  for 

\0i\<\,  0  <  &  <  1  +  20\ . 

At  the  stability  limit  0\  =  0i  =  2,  and  p  =  T~ 3. 


(3.216) 
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The  third  relaxation  scheme  considered  in  this  paper,  is  checkerboard  or  Red-Black  relaxation 
(RB).  It  consists  of  a  PJ  sweep  on  the  even  points  (k  =  0, 2, . . . ,  N  —  2),  an  update  of  the  solution 
and  residuals,  followed  by  a  similar  operation  on  the  odd  points  (k  =  1, 3, . . . ,  N  —  1): 


t>2i  =  (1  -  0)v2 k  +  0V2k-\ 

t>2*+l  =  (1  -  0)»2k+ 1  +  0V2 k 


k  =  0,1 . ±JV-1. 


In  the  Fourier  domain 


(3.22) 


1,  =  -  l0*f-1(l+f-1)vl+N/2,  l  =  0,l,...,N-l.  (3.23) 


This  implies  that  the  Fourier  components  at  /  and  l+N/2  are  coupled.  If  the  elements  (l,l  +  N/2) 
are  grouped,  then 


ft5'  )=Skf.6'  ),  Pm(iu  /  =  0,...  AN-l, 

V  Vl+N/ 2  )  V  v/+AT/2  /  \  *21  *22 )  2 


where 


and 


in(f)  =  l-0(l-f-1)-sn(-f), 

*i2(t)  =  -^2T-1(l  +  f-1)) 

*21  (T)  =  Sl2(-T), 

«22(T)  =  «n(-T).  • 


(3.24a) 

(3.246) 

(3.24c) 


Here  we  have  used  the  fact  that 


Ti+n/2  =  —Ti,  Tt  =  T  =  exp(2*il/N).  (3.25) 

Because  of  this,  the  expressions  for  s2i  and  *22  are  immediately  recovered  from  <u  and  *12-  The 
product  of  two  matrices  with  property  (3.24c)  again  has  this  property. 

The  eigenvalues  of  S  are 

=  |(in  +  S22)  d:  j[(«n  —  *22)*  +  4ii2«2i]^,  (3.26) 

Because  the  Vk  are  real,  we  have  |«i|  =  |6jv-j|,  so  only  the  values  l  =  0, 1, ... ,  or  6  £  [0,  t] 
have  to  be  considered.  Furthermore,  if  the  vi  are  grouped  according  to  (/,  l  +  N/ 2),  it  is  sufficient 
to  consider  /  =  0, 1, . . . ,  N/A  ox  6  £  [0,  x/2],  because  |v(+Ar/2|  =  i^Ar/2— il-  If  S  is  diagonal,  then, 
for  a  given  9  £  [0,  r/2],  the  modulus  of  A_  =  en(9)  describes  the  amplification  of  the  long  waves 
(0  £  [0,  x/2]),  whereas  the  modulus  of  A+  =  *22(0)  —  «n(»  —  9)  describes  the  growth  of  the 
short  waves.  If  S  is  not  diagonal,  the  eigenvalues  can  no  longer  be  assigned  to  these  parts  of  the 
spectrum. 

Following  [3],  the  best  relaxation  schemes  in  a  multigrid  context  are  those  for  which 

F  =  |p(0)|  (3.27) 

»€[**,*] 

is  as  small  as  possible.  Here  it  is  assumed  that  S  is  diagonal.  For  non-diagonal  S,  [3]  defines 

F  =  4  max  |s22(0)i  =  max  |sn(0')l-  (3-26) 

»6[0,*/2]  J'€[ir/2,ir] 

The  parameter  0,  or  the  pair  (/?i,/?2),  can  be  tuned  to  minimize  Jt.  Table  1  lists  the  optimum 
choices.  The  optimum  0  for  RB  is  the  root  of  0a  +  0  =  1.  However,  for  a  purely  convective 
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equation,  these  choices  may  not  lead  to  fastest  MG  convergence,  as  will  be  seen  in  the  next 
section.  The  main  point  is  the  diffusion  of  the  schemes  with  minimum  Ji.  It  will  be  shown  that 
smoothing  the  error  is  sufficient.  For  RB,  this  can  be  done  without  damping. 

We  continue  with  the  smoothing  rates  for  second-order  schemes.  The  differential  operator 
for  the  linear  case,  in  a  smooth  region  of  the  flow,  is 

Lh  =  i(l-T-1)(l+i(f-T’-1)).  (3.29) 


The  Point- Jacobi  scheme  has  a  growth- factor 


(3.30) 


There  is  no  value  /3  ^  0  for  which  this  scheme  is  stable. 

The  MS  scheme  with  second-order  spatial  accuracy  is  stable  for  0  <  0\  <  5, 0  <  02  <  2/?i . 
The  smallest  value  for  Ji  is  obtained  for  0i  =  |(1  +  3v^5),  02  =  (23  —  2y/S). 

Finally,  second-order  Red-Black  relaxation  is  stable  for  0  <  0  <  1.  The  smallest  value  of  Ji 
is  obtained  for  the  root  of  45/33  +  4/3  =  8.  A  variation  of  Red-Black  relaxation  is  obtained  if  the 
linear  distribution  characterized  by  Aj,  is  frozen  during  relaxation.  This  scheme,  however,  can 
not  be  stabilized  for  any  0  0  (contrary  to  the  remark  in  [16]). 

The  values  of  0  for  which  Ji  is  smallest  are  listed  is  Table  1. 

4.  Multigrid  relaxation 

4.1.  Coarse-grid  correction  operator.  The  convergence  speed  of  a  relaxation  scheme  can 
be  improved  by  the  use  of  coarser  grids.  The  communication  between  grids  is  controlled  by 
restriction  and  prolongation  operators.  Here  we  use  the  same  operators  as  in  [10].  Restriction 
brings  the  averages  of  pairs  of  cells  to  a  coarser  grid,  whereas  prolongation  brings  coarse-grid 
information  back  to  the  fine  grid.  In  this  section  we  will  analyze  the  behaviour  of  the  multigrid 
scheme  by  two-level  analysis  (cf.  [3]).  It  is  assumed  that  the  equations  on  the  next  coarser  grid 
(having  N/2  cells)  are  solved  exactly. 

The  two-level  amplification  operator  is 

M  =  (Shy>[I-  Lh](Shy\  (4.1) 

It  describes  the  result  of  applying  the  relaxation  operator  Sh  for  V\  steps,  restrict  the  residuals  on 
grid  h  to  the  coarse  grid  H  by  the  restriction  operator  l{f ,  solve  the  coarse-grid  problem  exactly, 
prolongate  the  coarse-grid  correction  to  the  fine  grid  by  /£,  add  this  correction  to  the  fine-grid 
solution,  and,  finally,  perform  1/2  additional  relaxation  sweeps.  We  will  now  derive  an  explicit 
expression  the  coarse-grid  correction  (CGC)  operator 


=  [/-/£(L") 


ULH)~ll"Lh] 


given  our  specific  problem. 

The  restriction  operator  I™  lets 

rlh  =  i(r2\  +  r$i+1),  k  =  0, 1, . . . ,  $N.  (4.3) 

Here  the  choice  H  =  2h  has  been  made.  Note  that  the  definition  of  7jJh  is  consistent  with  the 
definition  of  the  discretization  operator  Ih  in  (2.4),  and  preserves  the  conservation  form  of  the 
residual  (2.6).  The  prolongation  operator  is  ^2h  -  I:  the  coarse-grid  correction  €wis  simply 
distributed  over  the  fine  grid. 
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With  the  above  restriction  operator,  the  coarse-grid  residual  becomes 


r2kh  =  (llhLhvh)k  =  ±{T-T-')v2k,  k  =  0,1, . . . ,  ijV  —  1.  (4.4) 

In  the  Fourier  domain 

rfh  =  -  T"  *)(»,*  -  /  =  0, 1, . . —  1.  (4.5) 

The  result  is  0  for  the  shortest  and  longest  wave.  For  the  shortest  wave  T  =  —  1  and  the  restriction 
operator  makes  this  one  invisible  on  the  coarse  grid.  The  longest  wave  with  T  =  1  is  already 
made  invisible  on  the  fine  grid  by  the  residual  operator  Lh.  We  thus  have 

K(6)=l,  for  0  =  0,*.  (4.6) 

The  coarse  grid  operator  L2h  =  ^(1  —  T-2),  so  in  the  physical  domain  the  CGC  operator  K  lets 


V2i  =  «2 k  -  *>2i  +  l 

V2i+1  =  0 


Jfc  =  0,l,...,iJV-  1. 


(4.7) 


Here  the  longest  and  shortest  waves  have  been  ignored.  In  Fourier  space,  the  exact  solution  of 
the  coarse-grid  equation  yields 


-2j»  f  T(vf  /  =  1, . . . ,  1, 

'  \0,  /  =  0,  .  (4  8) 

The  prolongation  operator  ffh  lets 

ij*  =  t),h-|(l+f-1)02fc,  l  =  Q,l,...,N  -  1.  (4.9) 

Note  that  vfh  has  a  period  N/ 2,  whereas  if  has  a  period  N  in  /.  For  the  symbol  K  of  the  CGC 
operator,  the  elements  (1,1  +  N/ 2)  or  (6,6  +  *)  can  again  be  grouped  in  pairs,  as  in  (3.24a), 
yielding 

*n  =  i(l-r),  ku=l(l  +  f),  T  =  exp(»0),  0e(O,*/2].  (4.10a) 

The  elements  kji  and  kj i  are  found  through  property  (3.24c). 

Equation  (4.6)  states  that  the  coarse-grid  correction  has  no  effect  on  the  shortest  and  the 
longest  waves.  The  shortest  wave  can  be  taken  care  of  by  the  relaxation  operator  Sh .  The  longest 
wave  represents  a  constant  value.  If,  as  in  practical  computations,  a  sequence  of  coarser  grids  is 
used  instead  of  only  one,  this  constant  value  can  be  introduced  on  the  coarsest  grid  as  a  global 
constraint.  In  our  examples,  we  use  a  coarsest  grid  of  2  cells,  with  values  uj*  and  «**,  which 
should  obey  (2.14).  The  constant  C  =  c o  for  example  (2.2);  C  =  0  for  example  (2.3).  The  actual 
values  tij!  and  tij1  are  determined  as  the  exact  solution  on  the  grid  with  cell-size  hi,  which  obeys 
(2.13)  and  (2.14).  In  this  way,  we  effectively  have  K  =  0  for  9  =  0,  or,  more  precisely, 

*=(o  l)’  f°r  (*•  f  +  N/2)  =  (0, N/2)  or  (6,  *  —  0)  =  (0, *).  (4.106) 

Apart  from  this,  the  two-level  analysis  provides  a  reasonable  approximation  to  the  amplifica¬ 
tion  matrix  of  a  many-level  multigrid  correction,  because  the  use  of  multigrid  relaxation  for  the 
problem  on  the  grid  with  cell-size  2 h  usually  provides  a  very  good  solution. 


For  a  second-order  scheme  we  find  the  following  expression 

,  1 
'  1  + 


*11  =  i  -  4(i  + 


L  _  1/1  .  /Ti\  1  — i(^’  — ^  *) 

*12  -  j(l  +  T>  1+i(T*-T-i)’ 

f  =  exp(«0),  6  6  (0,  ir/2]. 


(4.11) 


For  6  =  0  we  again  have  (4.10b).  Alternatively,  we  may  consider  a  second-order-accurate  dis¬ 
cretization  on  the  fine  grid,  and  a  first-order  discretization  on  the  coarse  grid.  Then 


*n  =  i-i(i  +  T)  (l  +  Kr-r-1)) , 
*12=  5(1  +  r)  (1  —  \{t  —  r-1)j . 


(4.12) 


For  S  =  0we  should  use  (4.10b). 


4.2.  Asymptotic  convergence  rates.  The  quality  of  a  relaxation  scheme  in  combination  with 
the  CGC  operator  can  be  characterized  by 

A  =  max  p(M)-  (413) 

«€to,»/2r' 

Here  M  is  again  ordered  in  the  pairs  (6,  it  —  6),  so  that  the  maximum  is  effectively  taken  over  the 
entire  spectrum.  The  value  of  A  will  depend  only  on  v  —  v\  +  Vi ,  for  a  given  relaxation  scheme 
[13].  We  will  now  re-evaluate  the  three  relaxation  schemes  in  this  context. 

Consider  the  PJ  scheme.  The  smallest  value  of  A  is  obtained  for  0  =  |.  In  that  case  we 
have,  for  v  =  1  and  p—1, 

KhSh  =  -iCf-T-1)  (j  “}).  (4.14) 

The  spectral  radius  A=0.  Table  2a  shows  the  values  of  A  for  various  0.  For  a  second-order  scheme, 
PJ  is  unstable  for  long  waves,  but  the  coarse-grid  correction  stabilizes  these  waves.  This  scheme 
can  not  be  called  robust,  but  one  may  get  away  with  the  instability  in  practical  applications. 
For  a  first-order-accurate  discretization,  the  choice  0  =  \  and  u  —  1  is  clearly  the  best.  For 
second-order  accuracy,  the  first-order  CGC  yields  the  best  results.  It  does  not  pay  to  choose 
v  —  2,  because  the  improvement  of  the  damping  rarc  does  not  compensate  the  additional  work 
involved,  as  compared  to  v  =  1. 

The  results  for  the  Multi-Stage  scheme,  with  two  stages,  are  list’d  in  Table  2b.  For  a 
first-order-accurate  discretization,  one  better  uses  a  PJ  scheme.  For  second-order  accuracy,  a 
second-order  CGC  (p  =  2)  and  v  =  1  is  the  best  choice  in  terms  of  the  amount  of  work  involved. 

The  third  relaxation  scheme  to  be  considered  is  Red-Black.  Here  there  is  a  pitfall.  So  far, 
the  convection  speed  a  has  been  assumed  to  be  positive.  In  the  cases  considered  above,  the  final 
results  are  the  same  for  negative  a.  With  RB  relaxation,  the  asymmetry  in  the  relaxation  scheme 
interacts  with  the  asymmetry  in  the  CGC  operator,  thus  producing  different  results  for  positive 
and  negative  convection  speeds.  Instead  of  a  >  0  and  a  <  0,  we  can  just  as  well  consider  the 
effect  of  the  order  in  which  RB  is  carried  out.  So  far,  we  have  used  the  order  (even,  odd)  for  a 
grid  with  cells  *  =  0, 1, . . . ,  N  —  1.  The  final  results  for  a  <  0  and  the  order  (even,  odd)  are  the 
same  as  those  for  a  >  0  and  the  order  (odd,  even).  Table  2c  lists  the  results  for  both  the  order 
(even,  odd)  and  (odd,  even),  assuming  a  >  0.  Note  that  p(v  =  2)  ^  [p(v  =  l)]2,  in  contrast  to 
the  two  previous  relaxation  schemes. 
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The  first-order  scheme  with  0  —  1  produces  a  matrix  KhSh  =  0,  but  only  with  the  order 
(even,  odd)  and  a  >  0.  We  can  obtain  KhSh  =  0  for  positive  and  negative  convection  speeds,  if 
the  order  of  the  RB  scheme  is  adapted  to  the  direction  of  the  convection  on  the  coarser  grid.  If 
the  average  convection  speed  in  the  pair  of  zones  (2Jk,  2fc  +  l),  Jb  =  0,1,...,  ^N  —  l,  is  positive,  first 
the  even  and  then  the  odd  cells  are  relaxed,  whereas  for  a  negative  convection  speed  the  order 
(odd,  even)  is  used.  We  propose  to  name  this  scheme  Convective  Red-Black  (CRB)  relaxation. 
This  scheme  is  ideal  when  used  with  a  first-order-accurate  discretization,  because  the  combination 
with  the  CGC  operator  vanishes. 

4.3.  Ideal  relaxation  schemes.  RB  relaxation  is  known  to  be  order-free.  This  property  is 
lost  when  RB  is  embedded  in  our  multigrid  scheme.  CRB  effectively  is  a  marching  scheme  in 
the  multigrid  context,  like  Gauss-Seidel  (c.f.  [10,16]).  The  obvious  reason  for  this  behaviour  is 
the  structure  of  K  and  K  in  Eqs.  (4.7)  and  (4.10),  respectively.  Given  these  expressions,  the 
problem  of  designing  a  good  relaxation  scheme  can  be  approached  from  the  opposite  direction: 
which  scheme,  when  used  in  combination  with  (4.7)  or  (4.10),  gives  the  smallest  value  of  A  for  a 
given  amount  of  work? 

An  ideal  relaxation  scheme  obviously  should  have  the  property 

V2i  =  t»2t+i,  for  1  =  0,1,...,  \N  —  1.  (4.15) 

This  means  that  the  remaining  error  can  be  represented,  and  consequently  be  solved  exactly,  on 
the  coarse  grid.  CRB,  with  0  =  1  has  this  property.  A  more  economic  scheme  is  obtained  if 
relaxation  is  carried  on  half  the  number  of  cells.  Such  a  scheme  may  be  called  Red  or  Black,  but 
the  name  Semi  Red-Black  is  probably  less  confusing.  There  are  several  variants:  we  can  relax 
only  the  even  points,  or  only  the  odd,  or  take  the  direction  of  the  flow  into  account.  Assume 
a  >  0.  Then  the  order  (even,  odd)  results  in  a  relaxation  operator,  that,  for  0  =  1,  is  given  by 

&  1  +  f*1 

Sf  ~  2  ^  _i  +  t~1  l—T~1)'  (4-16) 

The  subscript  /  stands  for  forward  relaxation,  in  the  same  direction  as  the  flow.  The  backward 
variant,  corresponding  to  the  order  (even,  odd),  is 

q_l  (1  +  T-1  1+f-1)  ..... 

l-T-M’  (417) 

again  for  0=1.  The_ product  of  either  relaxation  operator  with  the  first-order  CGC  operator 
K  has  an  eigenvalue  A  =  0.  However,  for  one  combination  we  have  the  stronger  result  that  the 
product  vanishes,  namely: 

KSt  =  0.  (4.18a) 

The  effect  of  S»  is  equivalent  with  (4.15).  Furthermore, 

5;R  =  0  forf^l,  forf=l.  (4.186) 

If  the  longest  waves  (f  =  1)  are  ignored,  the  CGC  operator  K  causes  the  error  to  obey  (4.7), 
i.e.,  the  odd  cells  have  an  error  t>§*+i  =  Applying  Sj  on  the  even  cells  makes  the  error  in 
cell  2Jt  equal  to  the  error  in  cell  2Jb  —  1  by  convection,  resulting  in  a  zero  error  everywhere.  For 
v\  =  V2  =  1  we  have  SjKSj  =  0,  S/KSt  —  0,  and  StKSt  =  0.  However,  SiKSj  has  A  =  2  and 
is  useless. 


The  convective  variant  of  SRB  is  denoted  by  CSRB.  We  choose  the  combination  (4.18a)  if 
(yuvi)  =  (1,0),  and  (4.18b)  if  ( v\,vi )  =  (0, 1).  If  both  a  >  0  and  a  <  0  occur,  the  operators  Sj 
and  St  are  different  from  the  expressions  in  (4.16)  and  (4.17),  but  still  are  directed  along  with, 
or  against  the  flow.  For  instance,  Sj  on  the  pairs  (2k, 2k  +1),  £  =  0,1,...,  IjV  —  1,  updates  the 
even  cells  if  the  average  value  of  the  convection  speed  |(a*  ■+■  a*+1)  is  positive,  and  the  odd  if  this 
average  speed  is  negative.  The  operator  St  does  the  opposite. 

For  second-order  spatial  accuracy,  there  is  no  obvious  relaxation  scheme  that  removes  the 
error  completely  in  combination  with  the  CGC  step. 

It  should  be  noted  that  optimizing)!,  as  suggested  in  [3],  does  not  necessarily  give  the  optimal 
value  of  A. 

4.4.  Full  Multigrid.  The  performance  of  these  relaxation  schemes  has  been  tested  as  part  of  a 
Full  Multigrid  Scheme  (cf.  [3]).  We  start  on  a  grid  m  =  1  with  2  cells,  having  a  cell-size  hi,  were 
the  solution  is  computed  directly  from  the  nonlinear  discrete  equations,  using  the  conservation 
condition  (2.14)  as  an  additional  constraint.  On  the  finer  grids  m  =  2, 3, . . .,  with  N  =  2m,  we 
obtain  an  initial  guess  of  the  solution  by  interpolation  of  the  coarser  solution  (m—  1)  to  the  actual 
fine  grid  m.  The  corresponding  operator  will  be  denoted  by  Iff .  Then  a  number  of  nonlinear 
multigrid  or  FAS  (Full  Approximation  Storage)  cycles  are  carried  out  to  fiud  the  steady  state 
within  truncation  order. 

The  interpolation  of  the  final  solution  on  grid  m  —  1  to  grid  m  requires  an  order  of  accu¬ 
racy  consistent  with  the  solution.  We  first  describe  the  interpolation  of  a  second-order-accurate 
solution.  In  smooth  regions  of  the  flow,  the  solution  on  grid  m  —  1  is  assumed  to  be  piecewise 
parabolic: 

u?(z)  =  «?  +  *A<1>  +  -  &)A<2),  t  S  (4.19) 

Here  A^  equals  the  earlier  A*  of  Eq.  (2.10),  whereas  A^  is  a  discrete  approximation  to 

*i2>=  /  (4-20) 

The  interpolation  to  the  grid  with  cell-size  h  =  Hf  2  is  done  by 


u*  =  I^(nH) 


<4  =hj  dX  “ * ^  =  u?~  <  At1>’ 

*?+h 

«4+i  =  £  J  dxv?(x)  =  u?  +  LAll), 


Here  the  second-order  term  drops  out.  In  the  first-order  case  the  same  result  is  obtained.  In 
practice  A^  on  the  coarser  grid  with  cell-size  H  is  computed  by  means  of  the  same  averaging- 
limiting  procedure  as  before  (2.12). 

For  completeness  we  present  our  restriction  and  prolongation  operator  again.  Restriction 
involves  the  computation  of 

»k  =  j(«2t  +  «2i+l).  ' 

«f  =  £(4+4+i).  *  =  0,l,...,iiV-l.  (4.22) 

rh,k  =rk  ~  r(«t  )• 
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The  cosrse-grid  problem  now  becomes 
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rf  (tiH)  +  Tjffi  =  0,  or  Lf  (u")  =  sj?  +  r^t. 


(4.23) 


The  approximate  solution  uk  of  this  coarse-grid  problem  is  prolongated  back  to  the  fine  grid  by 


u2k 


-  .A 


=  “5*  +  («f  -  5(“2i  +  U2k+l))  =  5(«2*  -  u2t+l)  +  . 

*4+1  =  *4+1  +  («f  ~  5 (*4  +  *4+l))  =  5(U2i+l  "  «4)  + 


(4.24) 


In  the  actual  computer  program,  we  have  multiplied  s  J  and  rj  by  h,  which  simplifies  the  restriction 
of  the  residual  somewhat,  and  is  generally  recommended  for  non-uniform  grids  (cf.  [11]). 

The  multigrid  strategy  is  specified  by  i/\,  v?,  and  7.  Here  V\  denotes  the  number  of  sweeps 
before  restriction,  1/3  after  prolongation,  whereas  7  is  the  number  of  multigrid  cycles  on  the  coarser 
grid  H  before  returning  to  h.  A  flow-chart  for  this  strategy  can  be  found  in  [13,  pg.  45].  For 
7  =  1  we  have  a  V-cycle,  for  7  =  2  a  W-cycle.  A  less  expensive  version  of  the  latter  is  the  F-cycle 
[3],  also  with  7  =  2.  Here  the  first  cycle  from  grid  H  is  a  W-cycle,  the  second  a  V-cycle,  and  so 
on  for  all  coarser  grids. 


4.5.  Work.  The  amount  of  work  w‘  involved  in  one  multigrid  correction  cycle  on  grid  m  with 
2m  cells  is  given  by 

u>'  =  Am(7)  [uwx  +  wr  +  wp]  +  Bm(y)we.  (4.25a) 


Here  wx  is  the  amount  of  work  required  for  a  relaxation  sweep,  wr  for  restriction,  and  wp  for 
prolongation.  The  cost  of  the  exact  solution  on  grid  m  =  1  with  N  =  2  cells  is  expressed  by  we. 
The  constants  Am(y)  and  flm(7)  depend  on  the  total  number  of  grids  m  and  the  type  of  strategy: 


Amin)  =  E(7/2)n,  Bm(y)  =  \iy/2)m-\ 


(4.256) 


n=0 


Both  Am (7)  and  Bm( 7)  are  zero  for  m  <  2.  For  we  this  implies  that  the  cost  of  the  exact  inversion 
is  only  counted  during  the  multigrid  correction  cycle.  For  large  m,  we  have  Am(l)  ~  2,  yielding 
a  total  amount  of  work  proportional  to  N.  A  W-cycle  has  Am(2)  =  m  —  1  =  log2(AT/2).  For 
7  >  2  and  large  m,  we  have  Am( 7)  oc  Arl°8a7.  An  F-cycle  has  7  =  2,  but  still  results  in  an  O(N) 
algorithm.  It  has: 

^  =  E(n  +  14)n-  *£  =  ("»— lKj)"*-1.  (4.25c) 

n=0 


For  large  m,  A£  ~  4,  which  is  twice  the  cost  of  a  V-cycle. 

In  our  case,  the  evaluation  of  the  residual  is  the  most  time-consuming  part  of  the  algorithm. 
If  we  set  wc  =  1  for  the  PJ  scheme,  then  wt  =  2  for  MS.  For  RB  we  also  have  wz  =  2,  but  in 
other  cases  we  may  obtain  =  1  (for  instance,  if  first-order  flux- vector  splitting  for  a  system 
of  equations  is  used;  cf.  [9,10,11,16]).  Restriction  requires  the  computation  of  the  fine-to-coarse 
defect  correction,  which  involves  the  computation  of  the  residual  on  the  fine  and  the  next  coarser 
grid.  The  coarse-grid  residual  can  be  used  for  the  relaxation  scheme  or  further  restriction,  so  we 
only  have  to  count  the  fine-grid  work,  implying  wr  2;  1.  One  may  exploit  the  special  structure  of 
the  residual  to  lower  wT,  but  this  has  not  been  done  in  this  paper.  The  cost  of  prolongation  {wp) 
is  neglected.  The  first  time  the  exact  inversion  on  the  grid  with  N  =  2  cells  is  carried  out,  the 
discrete  source  term  can  be  used  directly  and  no  residual  has  to  be  computed.  In  this  case  wt 
is  neglected,  as  already  mentioned  above.  During  the  multigrid  correction  cycle,  we  need  rjfk  on 
the  coarsest  grid,  implying  that  the  residual  has  to  be  computed.  In  that  case  we  let  we  =  1. 
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The  total  amount  of  work  required  to  obtain  a  solution  on  grid  M  using  the  FMG-scheme  is 
given  by 

%  =  «[CW(7)(j/tu*  +«ir  +  wp)  +  Dm  (-y)«»«],  (4.26a) 

Here  «c  denotes  the  number  of  multigrid  correction  cycles  used  to  find  the  steady  state  on  each 
grid.  The  constants  Cm (7)  and  Dm (7)  depend  on  the  strategy  used: 

Af  M 

Cm(7)  =  J2  ^m(7)2m-",  DM(y)  =  Bm(7)2m~M.  (4.266) 

m= 2  m=2 

For  V-,  W-,  and  F-cycles,  this  results  in 

Cm(  1)  =  4  -  (i)w-J(Af  +  1),  DU{  1)  =  -  1); 

CM(2)=2(M-2)  +  (i)Af-2,  Z)Af(2)=l-(i)M-1;  (4.26c) 

Cm  =8-(i)"-1(M2  +  3M  +  4),  Z&  =  M(M  -  l)(i)M 

4.6.  Solution  error  and  stopping  criteria.  There  are  several  ways  to  determine  whether 
or  not  the  solution  on  a  given  finest  grid  has  converged.  A  simple  approach  is  to  compute  the 
number  of  cycles  required  for  convergence  below  the  discretization  error  from  the  estimated  A, 
and  hope  for  the  best.  In  this  case  the  initial  residual,  after  in  interpolation  of  the  stationary 
solution  on  the  coarser  grid,  should  be  brought  down  by  a  factor  (|)1+p- 

An  adaptive  approach  involves  a  comparison  between  the  residuals  and  an  estimate  of  the 
discretization  error.  The  latter  is  defined  by 

ri  =  Lhk  (Ik(u))  -  It  (!(«))  =  -rk(ut).  (4.27) 

Here  u  is  the  exact  stationary  solution  of  the  differential  equation,  and  tJ^  =  Jh(u).  The  operator 
Ih  is  given  by  Eq.(2.4).  Note  that  the  steady  solution  uk  of  the  discrete  problem,  obeying 
r* (uh )  =  0,  differs  form  u*  by  0(hp).  The  corresponding  error  in  the  stationary  solution 

Eh  =  (4.28) 

Following  [3],  the  discretization  error  in  the  residual  Th  can  be  estimated  by  assuming  rk(x)  ~ 
c(x)hp,  where  c(x)  is  independent  of  h.  Then 

r2hW  -  rHh  =  -  Th .  (4.29 a) 

From  this  expression  we  can  estimate 

r\x)  ~  (£)' J2\(r2*).  (4.296) 

The  fine-to-coarse  correction  r2A  is  computed  during  restriction  (4.22).  Although  Eq. (4.29a) 
should  contain  the  r2h  computed  from  steady  solutions  on  grid  h  and  2 h,  the  actual  value  during 
the  multigrid  cycling  can  be  used  as  an  estimate. 

If  one  does  not  want  to  use  the  non-steady  estimate  for  rfik,  rk  can  be  estimated  in  the 
following  way.  Assume  that  we  have  found  a  steady  solution  on  grid  m  —  1  with  cell-size  2 h. 
Then  compute  r*k  on  grid  m  -  2  by  (4.22).  From  this,  r4h  follows  through  (4.29a).  Next,  we 
apply  (4.29b)  twice  to  find  subsequently  rik  and  Tk. 


For  the  interpolation  (4.29b)  we  have  to  use  the  same  limiting-averaging  procedure  (2.11) 
as  for  u,  because  the  discretization  error  rk  can  be  discontinuous.  For  a  first-order  Godunov 
scheme,  applied  to  the  discontinuous  problem  (2.3),  rk  is  discontinuous  at  the  sonic  point  and  at 
the  shock.  Roe’s  scheme  smooths  the  discontinuity  at  the  sonic  point,  from  one  to  two  cells.  The 
second-order  scheme  has  a  smooth  rk  =  0(h 2)  almost  everywhere,  except  at  the  shock,  where  rk 
is  discontinuous  and  has  a  number  of  0(h)  spikes. 

Once  rk  has  been  estimated,  in  either  of  the  ways  sketched  above,  we  can  test  convergence 
by 

mi  <  e||r*||.  (4.30) 

For  smooth  problems,  the  L\~,  Z»j-,  or  Loo-norms  will  be  comparable.  A  local  comparison,  denoted 
by  Lioc,  will  generally  require  smaller  residuals.  A  natural  choice  is  e  =  1,  although  a  smaller 
value  may  be  necessary  if  r2k  is  computed  from  the  unsteady  solution. 

A  safe  way  to  determine  convergence  is:  test  (4.30)  for  the  Lt-norm,  with  t  =  1,  check  if 
the  Li-norm  of  the  residual  has  dropped  by  a  factor  (|)p+1 ,  and  test  if  the  Loo-norm  of  the 
residual  has  dropped  by  a  factor  The  motivation  for  this  triple  test  lies  in  the  properties  of  the 
discontinuous  solution.  The  Li-norm  is  fairly  insensitive  to  the  local  behaviour  at  the  shock  for 
a  first-order  scheme,  and  can  be  accurate  even  if  the  steady  state  has  not  been  fully  reached.  For 
a  second-order  scheme,  however,  the  local  0(h)-spikes  in  rk  may  increase  the  global  Li-norm, 
especially  if  the  solution  has  not  converged  yet.  In  that  case  testing  the  drop  in  the  residual  is 
a  safeguard.  The  Loo-norm  measures  the  convergence  near  the  shock.  After  grid-refinement  and 
local  relaxation,  the  initial  residual  around  the  shock  may  still  be  relatively  large,  so  considering 
only  the  drop  in  the  Loo-norm  of  the  residual  may  result  in  a  premature  stop.  The  combination 
of  the  three  tests  gives  a  satisfactory  performance  (§5.2). 

The  error  Ek  in  the  converged  solution  uk  can  be  estimated  in  a  similar  way  as  rk: 

Eh  =  (tyl2\(E3k),  E3k  =  TZ^Elk,  E3k  =  u3h-I3k(uk).  (4.31) 

A  first-order  Godunov  scheme  causes  Ek  to  be  discontinuous  at  the  sonic  point,  and  Roe’s  scheme 
again  smears  out  this  discontinuity,  but  only  slightly.  For  our  problem,  the  error  is  almost 
continuous  at  the  shock,  apart  form  an  0(h)  spike.  The  second-order  Godunov  or  Roe  scheme 
has  Ek  —  0(h2),  except  at  the  shock,  where  a  discontinuity  and  one  or  more  spikes  occur.  These 
spikes  are  of  0(h3)  for  the  present  problem,  which  has  a  practically  flat  solution  on  both  sides  of 
and  near  the  shock.  In  general  they  may  become  0(h)  in  the  worst  case. 


4.7.  Choice  of  the  time-step.  So  far,  we  have  not  specified  the  construction  of  a  fixed  /?. 
Because  the  value  of  0  under  consideration  may  correspond  to  the  stability  limit,  this  construction 
should  be  carried  out  with  care.  We  use  an  explicit  scheme  with 
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al  =  max 
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4>k  *  =  1»  4>kl* k  —  —  1. 


(4.32) 


If  the  expression  for  the  average  convection  speed  aj  happens  to  provide  zero,  we  simply  set  it  to 
h. 

The  above  choice  (4.32)  breaks  down  at  the  sonic  point  if  Godunov’s  scheme  is  used.  This 
can  be  seen  as  follows.  Consider  the  residual  in  cell  k.  Let  the  state  u*_i  <  0,  u*+i  >  0,  and 
u*_i  <  ti*  <  u*+j.  For  positive  u*,  the  first-order  Godunov  scheme  yields  r*  =  s*  -  For 
negative  ut  we  have  r*  =  s*  +  iuj.  If  the  residual  is  evaluated  at  a  coarser  grid  H ,  then  the 
fine-to-coarse  defect  correction  rjf  may  be  thought  to  be  absorbed  in  s*.  Equation  (4.32)  now 
yields  At  =  l/|«t|  if  0  =  I,  thus  providing  Newton’s  method  for  solving  nt  =  0.  However,  if  st 
and  it*  have  different  signs,  then  u*  has  to  change  its  sign  before  the  solution  u*  to  r*(u*)  =  0 


can  be  found  (assuming  that  the  conditions  ut_i  <  w*  <  u*+i,  ut-i  <  0,  u*+i  >  0,  remain  valid 
all  the  way  toward  the  steady  state).  In  that  case,  Newton’s  method  becomes  unreliable,  as  u* 
may  become  very  small,  so  that  r*  hardly  changes  and  At  becomes  extremely  large,  thus  causing 
extremely  large  changes  in  the  solution.  For  this  reason,  we  prefer  to  use  Roe’s  scheme. 

4.8.  The  singularity  at  the  shock.  In  the  discontinuous  case,  we  expect  two  kinds  of  problems. 
First,  the  interpolation  operator  I*h  (4.21),  that  computes  an  initial  guess  on  grid  h  from  the 
steady  solution  on  grid  2k,  will  cause  a  local  error  of  0(1)  at  the  shock.  In  smooth  regions  the 
error  will  be  0(hp).  Therefore,  something  special  has  to  be  done  at  the  shock  to  reduce  the 
0(l)-error  to  0(hp),  in  order  to  take  full  advantage  of  the  multigrid  scheme.  One  could  use  a 
more  sophisticated  kind  of  interpolation,  but  it  is  easier  to  use  local  relaxation,  which  we  have  to 
use  anyhow,  as  will  be  seen  below. 

The  second  problem  is  the  singularity  at  the  shock  in  the  discrete  equations.  Even  under 
ideal  circumstances,  the  asymptotic  convergence  rate  X  for  a  firsUorder  scheme  can  at  best  be  4, 
if  there  is  a  shock.  To  assert  this,  assume  that  we  have  two  grids,  a  fine  grid  m  with  cell-size  h, 
and  a  coarse  m  —  1  with  cell-size  2 h.  Also  assume  that  the  coarse-grid  problem  is  solved  exactly, 
and  that  the  residuals  are  very  small  everywhere,  even  at  the  shock.  In  that  case  the  errors  t>£ 
will  be  small,  except  at  the  shock,  where  v1^  can  still  be  large.  Without  loss  of  generality  we 
may  assume  the  index  M  of  the  cell  that  contains  the  singularity,  to  be  even.  The  coarse-grid 
error,  after  restriction,  is  V\M  =  HVM+VM+l)^hVM-  The  exact  solution  of  the  coarse  grid  is 
prolongated  back  to  the  fine  grid,  resulting  in 

w  1  ..h  M  X 

VM  -  VM  2 VM  ~  2vM’  (a 

-  „*  _  1.A  - 

VM+1  —  VM+1  2VM  —  2VJU- 

If  the  new  uj,  still  obeys 

1.  witl>  «Af~i  >  0,  Um+i  <  0,  (4.34) 

then  only  ujf+1  will  result  in  fairly  large  residuals  at  M  and  M  +  1.  An  ideal  relaxation  scheme 
such  as  CRB  with  0—1,  will  first  bring  v^f+1  down  to  practically  zero,  causing  the  residuals 
at  M  and  M  +  1  to  vanish.  The  error  v#  =  will  be  unaffected.  Thus  the  asymptotic 
convergence  rate  A  =  A  in  this  case.  In  general,  we  will  have  A  >  A. 

There  remains  an  open  question:  how  can  the  coarse-grid  problem  be  solved  under  the  above 
assumptions?  The  residual  on  the  coarse  grid  m  —  1,  after  restriction,  is  small,  so  how  can  v\ 
be  found?  Moreover,  the  coarse-grid  equations  will  be  plagued  by  the  same  singularity.  To  resolve 
this  issue,  we  have  to  consider  a  sequence  of  coarser  grids.  The  coarsest  grid  with  N  —  2  will 
notice  a  deviation  from  the  conservation  condition  (2.14),  even  if  the  residuals  are  zero,  and  make 
the  appropriate  correction.  This  correction  is  prolongated  to  the  next  fine  grid  (N  =  4).  If  many 
multigrid  cycles  on  N  =  4  are  made,  the  error  will  finally  vanish  on  this  grid,  with  a  damping 
rate  that  at  best  is  By  induction  we  conclude  that  in  this  way  the  exact  solution  on  grid  m—  1 
can  be  found,  as  long  as  f  is  sufficiently  large. 

If  a  relaxation  is  scheme  is  used  that  is  not  ideal,  and  if  smaller  values  for  j  are  used,  we  may 
end  up  with  a  situation  where  the  error  at  the  singularity  creates  large  residuals  around  the  shock 
after  the  coarsest  ( N  =  2)  grid  has  been  visited.  In  that  case,  we  may  not  have  convergence  at  all. 
We  therefore  need  a  special  subroutine  to  remove  the  residuals  around  the  the  shock,  and  at  the 
same  time  find  the  proper  solution  at  the  singularity.  In  view  of  (4.33),  we  propose  to  use  local 
relaxation  (cf.  [2])  by  a  Switched  Evolution/Relaxation  Method  (SER),  with  local  conservation 
imposed  as  an  extra  condition.  Let  the  set  C  denote  the  L  cells  where  this  local  relaxation  is 
carried  out.  Here  L  is  generally  much  smaller  than  the  total  number  of  cells  N.  Because  the 


information  for  conservation  is  obtained  form  the  coarser  grid,  C  is  chosen  to  cover  a  number  of 
coarse-grid  cells.  This  implies  that  L  is  even,  and  that  the  first  element  of  C  corresponds  to  an 
even  index  k  on  the  fine  grid. 

The  SER  scheme  is  given  by 


r„,  m,  n  €  C. 


(4.35a) 


For  small  time-steps,  this  scheme  is  practically  explicit,  hence  time- accurate,  whereas  for  large 
time-steps  Newton’s  method  is  obtained.  The  choice 

Kt  =  !/“«)•  (4-356) 

makes  (4.35a)  into  a  SER  scheme.  Similar  definitions  for  systems  of  equations  can  be  found  in 
[9,10,11,16].  For  At  —*  oo,  the  system  (4.35a)  becomes  singular  at  the  shock  ( drn/duM  =  0),  so 
the  corresponding  equation  is  replaced  by  the  requirement  of  local  conservation  (even  for  small 
time-steps): 

L 

Y  (A«u)m  =  0,  if  n  =  M,  m,  M  €  C.  (4.35c) 

m=l 

In  the  present  test  problem,  we  choose  L  —  4.  The  set  C  is  determined  by  finding  the  cell  ku 
where  the  residual  has  its  maximum  absolute  value.  For  L  =  4,  the  first  element  of  C  corresponds 
to  the  largest  index  k,  that  is  even  and  smaller  than  ku-  For  these  cells,  the  linear  system  (4.35)  is 
set  up  and  solved.  The  solution  is  updated,  and  the  new  residuals  in  C  and  the  two  neighbouring 
cells  is  computed.  This  process  is  repeated  until  a  maximum  number  of  steps  has  been  reached, 
or  until  condition  (4.30)  is  satisfied,  or  until  the  damping  rate  of  ||r*||  in  C  becomes  practically 
1.  For  the  latter,  we  use  the  Loo -norm.  The  entire  routine  may  be  repeated  a  number  of  times. 
Because  L  <  N,  in  general,  the  cost  of  applying  this  routine  is  low  compared  to  a  full  relaxation 
sweep.  This  routine  is  used  after  the  interpolation  of  the  coarse-grid  steady  state  to  the  fine  grid 
(4.21).  It  is  also  applied  once  directly  after  prolongation,  on  every  grid. 

4.9.  The  sonic  point.  The  linear  analysis  based  on  a  constant  convection  speed  a  is  not  valid 
at  the  sonic  point.  A  different  approach  is  required  to  find  the  damping  rates  around  this  point. 
We  will  consider  the  Point-Jacobi  scheme  for  a  first-order-accurate  solution  as  an  example. 

Let  the  sequence  of  cells  (2k,  2k  + 1, 2k  ■+  2, 2k  +  3)  contain  the  sonic  point.  Renumber  them 
as  (0, 1, 2, 3).  The  steady  state  U*  is  written  as  a*,  and  the  error  v*  =  a*  —  u*  is  assumed  to  be 
small  with  respect  to  a*.  A  local  expansion  of  at  can  be  made  around  ay. 

a»  =  aj  +  (k  -  2)ha'7  +  0(h3),  with  a7  >  0.  (4.36) 

On  the  coarse  grid,  we  denote  the  various  quantities  by  capital  letters  rather  than  using  super¬ 
scripts  h  and  H .  Let  Ao  =  ^(ao  +  aj)  <  0  and  Ai  =  %(a7  +  03)  >  0.  On  the  fine  grid  there  are  2 
cases,  with  a  different  position  of  the  sonic  point.  In  case  1,  we  let  a2  <  0, 03  >  0,  in  case  2,  we 
have  aj  <  0,  aj  >  0.  It  is  convenient  to  write 

case  1  :  aj  =  W2,  (4.37a) 

case  2 :  aj=  fafioj.  (4.376) 


Note  that  in  §3.2  we  neglected  terms  of  O(h).  This  can  not  be  done  at  the  sonic  point,  because 
a3  =  0(h). 
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The  linearized  residuals  can  be  found  from  the  fully  nonlinear  equations  by  neglecting  terms 
of  0(v2).  For  instance: 


ro  =  «o  -  j(«i  -  tio)  ’ 

*o  =  5(<*i  “  a?) 

«t  =  a*  -  t>t 


=>  tq  ~  aiVi  —  ao vq . 


(4.38) 


A  relaxation  step  with  Point-] acobi  lets  ti  =  v  —  /3rk/al,  where  aj  follows  from  (4.32).  For 
instance,  aj|  =  max(— ao,  —a i)  =  — oq.  In  this  way  we  find  in  case  1: 
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(4.39a) 


and  in  case  2: 


S<2)  =  /  -  0 


z1 

3a  — 2 
4-3a 

o 

1 

o  o 

-3a 

2+3a 

0 

(4.396) 


Here  S  works  on  the  vector  (v0,  t>i,  i>2,  vs)T .  The  exact  solution  of  the  coarse-grid  equations  yields, 
in  case  1: 

Vo  =  [(5  -  a)(4  +  a)t>o  -  2a(l  -  a)t>2  +  (2  -  a)(3  +  a)u3]/16,  ^ 

V\  =  [(1  —  a)(4  +  a)t>o  +  2a(3  +  a)t>2  +  (2  —  a)(7  +  a)v3]/16. 

In  the  other  case: 


Vo  =  [(4  —  3a)(5  +  3a)«o  +  3a(l  +  3a)vt  —  (2  —  3a)(l  +  3a)t>2  -f  3(1  —  a)(2  +  3a)t>3]/16, 

Vi  =  [(4  —  3a)(l  +  3a)vo  —  9a(l  —  a)vi  +  3(1  —  a)(2  —  3a)t>2  +  (7  —  3a)(2  +  3a)vs]/16. 

(4.406) 

From  this  solution  the  new  errors  on  the  fine  grid  can  be  found,  e.g.,  vq  =  vo  —  Vo-  In  this  way 
we  can  compute  the  CGC  operator  K.  The  spectral  radius  A  of  K.S  depends  on  (a).  For  P  —  \ 
and  a  €  [0, 1],  we  have 

case  1 :  0.318  <  I  <  0.359,  (4.41a) 

case  2  :  A  =  0.5 .  (4.416) 


Recall  that  we  found  A  =  0  in  the  case  of  constant  a.  Thus,  the  asymptotic  convergence  rate  will 
be  dominated  by  the  slower  convergence  at  the  sonic  point. 

For  the  other  relaxation  schemes  we  may  expect  a  similar  slow-down,  or  even  instability.  One 
could  consider  a  special  treatment  of  the  sonic  point,  e.g.,  by  using  additional  local  relaxation. 
In  practice,  only  convergence  down  to  the  truncation  error  is  desired.  The  local  relaxation  that 
is  applied  for  a  number  of  times  right  after  grid- refinement  will  not  only  take  care  of  the  region 
around  the  shock,  but  also  of  the  region  around  the  sonic  point,  since  the  truncation  error  is  large 
over  there,  and  therefore  also  the  initial  residuals.  This  may  be  sufficient  to  avoid  a  slow-down 
of  the  overall  convergence  during  the  first  few  cycles. 

4.10.  Summary.  We  have  analyzed  the  performance  of  a  few  relaxation  schemes  in  a  multigrid 
context.  Standard  multigrid  analysis  can  be  used  for  the  smooth  regions  of  the  flow.  Convection, 
which  is  the  physical  property  of  the  differential  equation  in  smooth  regions,  allows  for  the 


construction  of  ideal  relaxation  schemes.  Here  the  multigrid  technique  works  very  well.  On  the 
other  hand,  there  are  problems  at  positions  where  the  convection  speed  changes  sign. 

The  upwind  differencing  causes  the  loss  of  one  discrete  equation  in  the  cell  that  contains  the 
shock.  This  reflect  the  physical  property  of  shocks  to  absorb  incoming  signals,  but  at  the  same 
time  makes  the  steady-state  problem  ill-posed.  The  result  is  slower  convergence,  as  argued  in 
§4.8.  The  correct  steady  state  at  the  shock  can  be  recovered  by  using  local  conservation.  The 
required  information  can  be  obtained  from  the  coarser  grid.  The  local  relaxation  (4.35)  handles 
this  situation. 

At  the  shock  we  lose  a  discrete  equation.  At  the  sonic  point,  we  effectively  gain  one,  as  is 
obvious  from  (4.39).  This  causes  the  local  damping  rate  to  deviate  from  its  overall  value.  If 
the  local  damping  is  larger  than  the  overall,  one  may  consider  a  special  treatment  of  the  sonic 
point.  Because  there  is  an  additional  discrete  equation,  one  could  find  the  local  steady  state 
directly  by  considering  two  cells  simultaneously  (but  only  in  the  first-order  case).  This  implies, 
again,  local  relaxation.  Alternatively,  the  relaxation  parameter  f3  could  be  optimized  separately 
for  the  sonic  point.  However,  we  have  decided  to  do  nothing  special  at  the  sonic  point  during  the 
multigrid  correction  cycle,  as  we  expect  that  the  local  relaxation  applied  after  grid-refinement 
will  be  sufficient  to  find  the  steady  state  within  truncation  error. 

5.  Numerical  examples 

5.1.  Asymptotic  convergence  rates.  The  above  algorithm  has  been  coded  for  the  smooth 
and  discontinuous  examples  of  Eqs.(2.2)  and  (2.3),  respectively.  For  the  smooth  case,  we  expect 
a  good  agreement  with  the  estimated  damping  rates.  In  the  discontinuous  case,  local  relaxation 
should  take  care  of  the  singularity  near  the  shock.  The  sonic  point  may  still  cause  deviations 
from  the  predicted  damping  rates  for  constant  convection  speed. 

The  asymptotic  convergence  rate  A  holds  for  a  large  number  of  points  N  and  an  exact  CGC- 
operator.  For  finite  N,  there  will  be  an  O(h)  deviation.  To  make  this  small,  we  have  made  runs 
for  a  grid  with  N  =  512  cells.  For  the  first-order  CGC,  we  use  the  CRB  scheme,  with  f)  sr  1, 
v\  —  0,  1/2  =  1,  and  7  =  2.  This  covers  the  first-order-accurate  solutions  [p  =  1],  and  the  second- 
order-accurate  ones  that  are  obtained  with  a  first-order  CGC  [p  =  (2,  l)j.  For  the  second-order 
CGC,  we  use  the  same  relaxation  scheme  on  the  grid  m—  1  as  on  the  finest  grid  m,  with  the  same 
i/i  and  i/2  as  on  the  finest  grid.  On  grid  m  —  1,  we  let  7  =  3  to  make  sure  that  the  coarse-grid 
problem  is  solved  almost  exactly.  For  the  subsequent  coarser  grids  m  —  2,  m  —  3, . . .,  the  amount 
of  work  would  become  rather  large  with  7  =  3,  so  we  use  the  same  first-order  CGC  as  above,  but 
now  with  respect  to  grid  m  —  1  rather  than  to  the  finest  grid  m. 

Local  relaxation  is  only  applied  in  the  discontinuous  case.  After  grid-refinement,  at  most 
8  sweeps  are  made  for  the  first-order-accurate  solution,  and  at  most  12  for  the  second-order- 
accurate  one.  In  both  cases,  the  search  for  the  largest  residual  is  carried  out  at  most  4  times, 
so  for  the  first-order-accurate  solution  the  local  relaxation  is  repeated  at  the  same  4  points  for 
at  most  2  times,  for  the  second-order-accurate  solution  at  most  3  times.  During  the  multigrid 
correction  cycle,  we  only  carry  out  local  relaxation  once  on  4  cells. 

The  evaluation  of  the  numerical  damping  rate  A  involves  some  arbitrariness.  The  initial 
values  of  the  state  quantities  u*  are  taken  from  the  solution  on  the  coarser  grid.  Thus,  the  higher 
frequencies  will  dominate  the  error.  During  the  first  few  cycles,  the  damping  rate  is  mainly 
determined  by  the  behaviour  of  these  high  frequencies.  After  this  phase,  the  convergence  rate 
becomes  more  or  less  constant,  until  round-off  errors  start  to  dominate.  The  numerical  damping 
rates  have  to  be  estimated  from  the  part  of  the  convergence  history  between  the  initial  phase  and 
the  noise- level. 

For  the  results  shown  in  Table  3,  10  multigrid  cycles  on  the  finest  grid  were  carried  out,  from 
which  the  damping  rate  per  cycle  was  estimated.  For  the  larger  values  of  A,  this  usually  involved 
all  but  the  first  cycle. 


It  turns  out  that  the  actual  result  is  practically  the  same  for  the  L\-,  Li-,  and  Loo -norm  of 
either  the  error  vh  or  the  residual  rh,  although  the  Loo-norm  is  less  constant  from  step  to  step. 
The  results  in  Table  3  are  based  on  the  Li-norm,  and  are  given  for  both  the  error  vh  and  the 
residual  rh. 

The  linear  analysis  predicts  identical  results  for  (iq , 1/2)  =  (1,0)  and  (1/ 1,1/3)  =  (0, 1).  The 
differences  that  occur  in  Table  3  give  an  indication  of  the  accuracy  of  the  estimated  damping 
rates,  as  do  the  differences  between  the  results  for  vh  and  r*. 

In  most  of  the  cases,  there  is  a  good  agreement  between  analysis  and  experiment.  The  sonic 
point,  occurring  in  the  discontinuous  case  causes  slower  convergence  for  PJ,  p  =  1,  for  MS  and 
RB,  if  p  =  1  and  1/1  =  i/2  =  1,  and  for  SRB  and  CSRB  if  p  —  1.  The  SRB  and  CSRB  schemes 
are  mildly  unstable  at  the  sonic  point  for  second-order-accurate  solutions  [p  =  (2, 1)  or  p  =  2], 
although  not  in  all  cases.  This  instability  starts  to  dominate  after  the  rest  of  the  solution  has 
converged  toward  the  steady  state. 

5.2.  Fast  convergence.  The  above  results  show  that  our  linear  analysis  is  valid.  It  can  be  used 
to  evaluate  the  quality  of  the  several  relaxation  schemes.  For  practical  purposes,  these  asymptotic 
convergence  rates  are  less  important.  In  a  FMG-context,  the  initial  guess  is  obtained  from  the 
steady  solution  on  the  coarser  grid.  In  that  case,  the  damping  of  the  first  multigrid  cycle  can  be 
much  better  than  predicted  by  the  asymptotic  damping  rate.  The  initial  error  must  be  damped 
by  a  factor  (5)1+p  to  become  comparable  to  the  truncation  error.  For  some  relaxation  schemes, 
this  is  already  accomplished  before  the  asymptotic  regime  is  entered. 

Consider  the  computation  of  a  first-order-accurate  steady  state.  For  most  of  the  examples 
presented  in  Table  3,  the  asymptotic  damping  rates  A  <  0.25,  so  one  F-cyde  per  grid  can  be 
expected  to  suffice.  Table  4a  shows  results  obtained  in  this  way.  We  use  the  optimal  values  for 
the  relaxation  parameter  0,  as  in  Table  3,  but,  contrary  to  $5.1,  the  same  relaxation  scheme  is 
used  on  every  grid.  The  quantity 

^  =  ||«‘-uk||f/||«*-uj||f.  (5.1) 

Here  the  numerical  steady  state  u*  has  a  zero  residual,  whereas  u,  is  the  cell-averaged  stationary 
solution  of  the  differential  equation.  The  total  amount  of  work  w  is  computed  by  counting  the 
number  of  cells  in  which  the  residual  is  evaluated,  and  dividing  by  the  number  of  cells  on  the 
finest  grid.  For  the  smooth  solution,  these  values  are  consistent  with  (4.26).  The  total  amount  of 
work  w  for  the  discontinuous  case  is  higher,  because  of  the  local  relaxation  sweep  on  4  cells  after 
prolongation.  Because  F-cycles  involve  many  visits  to  the  coarser  grids,  the  relative  amount  of 
work  spent  on  local  relaxation  is  somewhat  large. 

Table  4a  shows  that  cycles  with  (1/1,  1/3)  =  (0,1)  are  preferable  above  those  with  (('1,1/3)  = 
(1,0).  The  smaller  amount  of  work  for  the  former  in  the  discontinuous  case  is  due  to  the  fact 
that  the  local  relaxation  routine  is  applied  in  the  neighbourhood  of  the  residual  of  largest  size, 
hence  requires  the  residual  to  be  known  everywhere.  When  uj  =  0,  the  residual  is  computed  on 
the  entire  current  grid,  although  it  is  modified  in  only  4  +  2  cells. 

The  Point-Jacobi  scheme  with  (iq,  1/3)  =  (0, 1)  is  clearly  the  most  cost-effective  choice.  The 
Semi  Red-Black  scheme  is  slightly  faster,  but  obviously  less  robust.  For  the  latter,  one  cycle  is 
not  always  sufficient. 

For  a  second-order-accurate  solution  the  use  of  a  first-order  CGC  appears  more  attractive 
than  a  second-order  CGC,  because  the  second-order-accurate  residual  is  substantially  more  ex¬ 
pensive  than  a  first-order-accurate.  Table  4b  shows  results  obtained  with  a  first-order  CGC,  using 
a  fixed  number  of  F-cycles  on  every  currently  finest  grid.  The  same  relaxation  scheme  is  used  on 
the  fine  and  coarser  grids,  but  with  a  different  optimal  parameter  0  for  the  finest  and  tjie  coarser 
grids.  The  number  of  cycles  k  is  based  on  the  estimated  asymptotic  convergence  rate  A.  For  PJ, 
MS,  and  RB,  we  obtain  convergence  well  below  the  truncation  error.  The  SRB  scheme,  which 


is  unstable  at  the  sonic  point  in  some  cases  ($5.1),  performs  less  satisfactorily.  In  some  cases  we 
need  more  cycles  to  obtain  a  converged  solution,  in  others  also  W-cycles  instead  of  F-cycles. 

5.3.  Deferred  Correction  technique. 

In  the  previous  paragraph,  a  first-order  CGC  was  used  to  obtain  second-order-accurate  solu¬ 
tions.  Instead  of  using  first-order-accurate  residuals  only  on  the  coarser  grids,  one  may  use  them 
on  the  finest  grid  as  well.  In  that  case  we  compute  a  term 

*k  =  *-£(p  =  2)  -  rj(p  =  1)  (5.2) 

at  the  begin  of  every  multigrid  correction  cycle.  The  term  is  added  to  the  first-order  residual 
during  relaxation  and  restriction.  This  technique  is  known  as  the  deferred  correction  technique 
(cf.  [3]).  The  computation  of  (5-2)  involves  some  additional  work,  but  this  is  compensated  by 
the  lower  cost  of  the  first-order  scheme. 

For  the  local  relaxation  after  grid-refinement,  sis  applied  in  the  discontinuous  case,  we  use  the 
second-order-accurate  residual.  The  quantities  <rj  are  computed  at  the  begin  of  every  multigrid 
correction  cycle,  and  kept  fixed  until  the  next  cycle.  Because  it  is  not  clear  what  the  asymptotic 
convergence  rate  will  be,  we  have  used  the  stopping  criterion  (4.30)  with  e  =  1,  and  with  the 
additional  requirement  that  the  Li-norm  of  the  residual  must  drop  by  a  factor  and  its  Loo -norm 
by  i. 

Numerical  results  are  shown  in  Table  4c.  Again,  PJ  with  (»'i,i'2)  =  (0, 1)  seems  to  be  the 
best  choice.  The  SRB  scheme  does  not  convergence  is  some  cases,  as  might  have  been  expected 
by  now. 

5.4.  Tau-extrapolation.  In  $4  it  was  discussed  how  to  estimate  the  truncation  error  rh.  These 
estimates  can  actually  be  used  to  improve  the  accuracy  of  the  solution.  The  technique  is  call 
r-extrapolation. 

The  simplest  implementation  estimates  r7h  from  rjjh  according  to  (4.29a),  and  adds  this  in¬ 
stead  of  T)Jh  to  the  coarse-grid  residual  (4.23).  However,  this  introduces  an  inconsistency  between 
the  finest  and  coarser  grids,  which  slows  down  convergence  when  the  residual  becomes  comparable 
to  the  truncation  error.  Consistency  can  be  recovered  by  adding  the  current  estimate  of  rh  to 
the  fine-grid  residual  during  relaxation.  The  initial  r*  at  the  beginning  of  the  iterations  on  grid 
h,  can  be  estimated  from  r^J,  using  (4.29).  At  the  restriction  from  h  to  2h  we  compute  r£h  and 
modify  the  residual  at  grid  2 h  by  means  of  r3\  At  the  same  time  a  new  estimate  of  the  fine-grid 
truncation  error  Th  is  obtained  through  (4.29b). 

Table  5  shows  the  errors  in  the  steady  solutions  obtained  in  this  way.  The  first  and  second 
row  give  the  results  without  r-extrapolation.  Both  the  L\-  and  Loo-norm  of  the  error  are  given. 
These  values  are  obtained  by  averaging  \\Eh\\/hp  from  the  grids  with  N  =  32  to  N  =  512.  The 
third  and  fourth  row  show  the  effect  of  r-extrapolation:  the  order  of  the  solution  is  increased  by 
1.  An  exception  is  the  discontinuous  case.  As  remarked  in  $4,  rh  is  discontinuous  at  the  shock 
and  at  the  sonic  point  if  p  =  1,  and  has  0(A)-spikes  around  the  shock  if  p  =  2.  Consequently,  the 
r-extrapolation  loses  its  accuracy  at  these  isolated  spots.  This  shows  up  in  the  Loo-norms  for  the 
discontinuous  case. 

6.  Conclusions 

The  analysis  and  experiments  show  that  the  standard  multigrid  techniques  can  be  applied  to  a 
nonlinear  hyperbolic  conservation  law  of  the  type  studied  here.  In  regions  where  the  convection 
speed  is  well  above,  or  below,  zero,  ideal  damping  rates  [of  0(h)]  can  be  obtained.  If  the  con¬ 
vection  speed  changes  its  sign,  however,  this  ideal  behaviour  is  lost.  At  the  shock  we  have  a 
singularity  in  the  discrete  equations.  This  is  a  result  of  dropping  the  time-accuracy.  The  sin¬ 
gularity  can  be  removed  by  augmenting  the  equations  with  Rankine-Hugoniot  jump  conditions, 


which  is  equivalent  to  requiring  local  conservation  as  an  additional  constraint.  Information  about 
conservation  is  obtained  from  coarser  grids,  through  a  special  kind  of  local  relaxation.  At  the 
sonic  point,  the  analysis  is  not  applicable,  but  in  most  cases  this  is  not  a  serious  problem. 

We  have  analyzed  3  schemes:  Point-Jacobi,  a  Two-Stage  scheme,  and  Red-Black  relaxation. 
Point- Jacobi  appears  the  simplest  and  most  cost-effective,  both  for  first-order  and  second-order 
accuracy.  Its  cheaper  variant,  Semi  Red-Black,  is  theoretically  appealing,  but  not  very  robust. 
Normal  Red-Black  relaxation  loses  its  direction-free  property  when  used  in  combination  with 
the  CGC  operator.  One  either  needs  underrelaxation  (0  =  5),  or  must  take  the  direction  of 
convection  into  account  (Convective  Red-Black). 

The  main  mechanisms  for  smoothing  are  damping  and  convection.  Damping  causes  smearing 
and  dissipation  of  the  error.  Convection  shifts  errors  around.  It  can  profitably  be  used  to  make 
errors  equal  in  neighbouring  cells,  which  is  ideal  for  a  smoother.  In  the  presence  of  a  shock  or 
a  boundary,  there  is  another  mechanism  to  get  rid  of  errors.  In  the  neighbourhood  of  a  shock, 
convection  shifts  errors  toward  the  shock,  where  they  are  absorbed,  whereas  near  a  boundary  the 
simply  leave  the  computational  domain.  Thus,  a  boundary  or  shock  will  change  the  local  the 
behaviour  of  the  multigrid  scheme.  It  may  be  that,  on  a  relatively  coarse  grid,  the  schemes  with 
underrelaxation  (0  <  1)  are  less  efficient  than  those  with  0=1,  even  if  the  latter  are  unstable. 
However,  on  much  finer  grids,  this  can  no  longer  be  true. 

First-order-accurate  steady  states  can  be  obtained  in  one  F-cycle,  independent  of  the  number 
of  cells  in  the  computational  domain.  A  second-order-accurate  solution  requires  about  4  such 
cycles,  using  either  a  first-order  CGC  scheme,  or  the  deferred  correction  technique. 

First-order  upwind  differencing  is  a  good  starting  point  for  the  construction  of  efficient 
solvers.  It  directly  reflects  the  physical  behaviour  of  the  differential  equation  and  does  not  suffer 
from  the  spurious  oscillatory  phenomena  that  require  the  use  of  artificial  viscosity  in  central- 
difference  schemes.  It  is  therefore  easier  to  pinpoint,  and  interpret,  the  difficulties  that  may 
occur  in  actual  applications.  In  addition,  efficient  relaxation  methods  are  easier  to  construct  for 
upwind-difference  than  for  central-difference  schemes. 

The  present  results  can  be  easily  generalized  to  systems  of  equations  in  one  dimension, 
specifically  the  Euler  equations  of  gas-dynamics.  Roe’s  scheme  allows  the  use  of  characteristic 
equations  without  violating  conservation.  For  each  of  these  equations,  the  above  results  are  valid. 
Even  the  convective  variants  can  be  generalized,  as  can  the  ”time-step”  (4.32). 

It  should  be  remarked  that  in  one  dimension  we  may  just  as  well  use  tridiagonal  solvers 
[9].  For  two  dimensional  problems  the  use  of  the  multigrid  technique  is  more  important.  Un¬ 
fortunately,  analyzing  relaxation  schemes  and  the  CGC  operator  is  much  harder,  because  the 
two-dimensional  Euler  equations  of  gas-dynamics  are  not  simultaneously  diagonalizable  in  both 
spatial  dimensions.  A  scalar  model  equation  of  the  form 

+  aui  +  feu,  =  0  (6.1) 

is  not  applicable  (except  when  both  the  x-  and  y-components  of  the  velocity  are  supersonic)  and 
one  needs  to  resort  to  systems  of  equations,  which  complicates  the  analysis. 
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Table  1.  Choices  of  0  for  which  p  is  smallest.  The  results  are  given  for  one  relaxation  step 
(v  —  1).  Second-order  PJ  is  unstable. 
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Table  2a.  Asymptotic  multigrid  convergence  rates  for  Point-Jacobi  relaxation.  Optimal  values 
are  denoted  by  an  asterisk.  Second-order  PJ  is  unstable.  The  notation  p  =  2, 1  means:  p  =  2  on 
the  finest  grid,  p  =  1  on  the  coarse. 
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'able  2b.  Asymptotic  two-grid  convergence  rates  for  Two-Stage  relaxation.  For  /?i  =  0  this 
:heme  reduces  to  Point-Jacobi. 
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Table  2c.  Asymptotic  multigrid  convergence  rates  for  Red-Black  relaxation.  The  value  of  A 
depends  on  the  order  in  which  the  relaxation  is  carried  out:  (e,  o)  means  that  first  the  even  and 
then  the  odd  points  are  relaxed  (k  =  0, 1, . . . ,  N  —  1).  For  (o,_e)  this  order  is  reversed.  All  values 
are  computed  for  a  >  O.If  a  <  0  then  the  values  A (e,o)  and  A(o,e)  should  be  interchanged.  The 
above  optimal  values  of  A  for  the  ordering  (o,  e)  are  also  optimal  if  both  a  >  0  and  a  <  0  occur. 
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Table  3a.  Asymptotic  multigrid  convergence  rates  for  Point-Jacobi  relaxation,  as  estimated  from 
the  numerical  experiments  on  a  grid  with  512  cells. 
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Table  3c.  Asymptotic  multigrid  convergence  rates  for  Red-Black  relaxation. 
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Table  4a.  Convergence  levels  and  amount  of  work  for  a  first-order- accurate  steady  state  on  a 
grid  with  512  cells,  when  one  F-cycle  is  carried  out  on  every  grid  (k  =  1,  7  =  2). 
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Table  4b.  Convergence  levels  and  amount  of  work  for  a  second-order-accurate  steady  state  on 
a  grid  with  512  cells,  when  k  F-cycles  are  carried  out  on  every  grid.  A  first-order  CGC  is  used 
(P  =  2,1). 
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Table  4c.  Convergence  levels  and  a mount  of  work  for  a  second-order-accurate  steady  state  on  a 
grid  with  512  cells,  when  the  deferred  correction  technique  is  used. 


